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We present a pedagogical account of the Lee- Yang theory of equilibrium phase transitions 
and review recent advances in applying this theory to nonequihbrium systems. Through 
both general considerations and explicit studies of specific models, we show that the Lec- 
Yang approach can be used to locate and classify phase transitions in nonequihbrium steady 
states. 

I. INTRODUCTION 

In this work we seek a mathematical understanding of phase transitions in the steady state of 
stochastic many-body systems. Systems at equilibrium with their environment provide examples 
of such steady states, and the mechanisms underlying equilibrium phase transitions are long known 
and understood. Experimentally, one can distinguish between two types of transition: the first- 
order transition at which there is phase coexistence, e.g. between a high-density solid and a low- 
density fluid, and the continuous transition at which fluctuations and correlations grow to such an 
extent as to be macroscopically observable. 

From a thermodynamical perspective, one can understand first-order transitions by associating 
with each phase a free energy. For a given set of external parameters, the phase 'chosen' by 
the system is that with the lowest free energy, and so a phase transition occurs when the free 
energies of two (or more) phases are equal. The sudden changes in macroscopically measurable 
quantities that take place at first-order transitions are described mathematically as discontinuities 
in the first derivative of the free energy. Discontinuities in higher derivatives relate to continuous 
(higher-order) phase transitions. 

The tools of equilibrium statistical physics allow one — in principle at least — to express the free 
energy solely in terms of microscopic interactions. More specifically, the free energy is given by the 
logarithm of a partition function, a quantity that normalises the steady-state probability distri- 
bution of microscopic configurations. Initially it was not universally accepted that this approach 
could faithfully describe phase transitions, in particular the first-order solid-fluid transition [1]. In 
order to show that the statistical mechanical approach can reproduce the correct discontinuities 
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in the free energy at a first-order transition, Lee and Yang [2] introduced a description of phase 
transitions concerning zeros of the partition function when generahsed to the complex plane of 
an intensive thermodynamic quantity. Initially, the theory was couched in terms of zeros in the 
complex fugacity plane which is appropriate for fluids in the grand canonical (fluctuating particle 
number) ensemble. However, by mapping a lattice gas onto the Ising model [3], the theory was 
found to hold equally well for a magnet in a fixed magnetic field. Later [4-6] it became clear that 
the distribution of zeros in the complex temperature plane can reveal information about phase 
transitions in the canonical ensemble. 

In section II we present a brief, self-contained discussion of the Lee- Yang theory of equilibrium 
phase transitions which relates nonanalytic behaviour in the free energy to zeros of the partition 
function. In the remainder of this article, we move our attention to phase transitions in nonequi- 
librium steady states, that is, those that carry currents of mass, energy or some other quantity. 
We give therefore a very brief introduction to the established practice of modelling noncquilibrium 
physics through stochastic dynamics in section III. In particular, we will identify a candidate quan- 
tity for use in place of an equilibrium partition function when performing an Lee- Yang analysis of 
phase transitions. Finally, in section IV we review recent work in which the Lee- Yang theory has 
been successfully applied to specific nonequilibrium phase transitions. 

II. OVERVIEW OF LEE- YANG THEORY OF EQUILIBRIUM PHASE TRANSITIONS 

In this section we recount the relationship of zeros of a partition function to phase transitions 
in the system it describes. We present here a brief but self-contained tour of the key points. For 
more detail (and mathematical rigour) the reader is referred to textbooks, such as [7] or, for a more 
advanced presentation, [8], as well as early treatments of the subject [2-6]. 

For concreteness we shall describe the theory in terms of a simple model system of N spins in 
thermal contact with a heat reservoir. In this model system the energy of the system can take 
values E = ne with n = 0, 1,2, ... ,M and the number of microstates corresponding to the n^^ 
energy level is g{n). The canonical partition function is given by 

M 

ZN{f3) = J29i^')^M-Pne) (1) 

n=0 

where (3 is the inverse temperature. To investigate the zeros of this partition function, it is useful 
to make the change of variable z = e"^^. Then, the partition function is explicitly a polynomial in 
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z of degree M and can be factorised as 
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Clearly the quantities Zn are the M zeros of the partition function Zn(z); meanwhile, k is a 
constant which we can safely ignore in the following. 

Since all the g{n) are positive, no zero of Zn{z) can be real and positive: that is, the zeros Zn 
will generally lie in the complex plane away from the physical values of z which lie on the positive 
real axis. We now define for all complex z except the points z = Zn the (complex) free energy 



and note that a Taylor series expansion of hj^{z) around a point z ^ Zn has a finite radius of 
convergence given by the distance to the nearest zero from z. This then implies that that hN{z) 
can be differentiated infinitely many times over any region of the complex plane that is devoid of 
partition function zeros. Since we identify a phase transition through a discontinuity in a derivative 
of the free energy, we see that such a transition can only occur at a point zq in the complex plane if 
there is at least one zero of the partition function Zn(z) within any arbitrarily small region around 
the point zq. 

Clearly this scenario is impossible if the number of zeros M is finite, except at the isolated 
points Zfi where the free energy exhibits a logarithmic singularity. Since such a point cannot lie 
on the positive real z axis, there is no scope for a phase transition in a finite spin system, such as 
the simple example (1). On the other hand, if the partition function zeros accumulate towards a 
point zq on the real axis as we increase the number of spins N to infinity there is the possibility of 
a phase transition. 

In order to deal with the thermodynamic limit (see [8] for rigorous considerations) we shall 
assume that the limit 



hN{z) = 



In Zn{z) 
N 



(3) 



Using (2) we rewrite this as 




(4) 



h(z) = lim 
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exists and we may write 




(6) 
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where p{z) is the local density of zeros in the complex- 2; plane in the thermodynamic limit. 

Since the imaginary part of this complex free energy h{z) is multi-valued it will at times be 
more convenient to work with the potential ip{z) defined as 



<p{z) = Rehiz) = J dz'p{z') In 1 - ^ . (7) 



An expression for the density p{z) in terms of the potential ip{z) is easily obtained once one 
recognises that In is the Green function for the two-dimensional Laplacian. Specifically, if 
z = X + iy, we have 

^ + ^) + = 27r6{x)5{y) . (8) 

Using this expression, we can take the Laplacian of both sides of (7) to find 

Pi^) = ^VV(^) ■ (9) 

Such an equation is familiar from electrostatics and so ip{z) is analogous to an electrostatic poten- 
tial. 

In analogy to electrostatics, as long as we can integrate p{x) over bounded regions containing 
any singularities of p, the potential will be a continuous function [9]. The significance of this 
statement is that we can derive a rule for locating phase boundaries given a partition function. Let 
us suppose that around points zi and Z2 in the complex plane, one has (analytic) expressions for 
the potential ipi{z) and (p2iz) such that (pi{z) ^ (p2iz) (i.e. not the same function over the entirety 
of the complex plane). In order for the potential to be continuous at all points on the complex 
plane, we must have a phase boundary at those values of z for which the condition 

Mz) = M^) (10) 

holds. This is basically the definition of a transition mentioned above in the introduction. Since 
(fi and ip2 are different functions, some derivatives of ip{z) will not exist at these values of z and 
we expect the density of zeros at these points to be non-zero. (It is also possible for zeros to be 
present at other points in the complex plane, but we do not need to consider this possibility here.) 

Typically, a solution of (10) describes a curve C that intersects the positive real z axis at a 
point zq. Having already established that we require zeros to accumulate at the point zq in the 
thermodynamic limit for a physical phase transition to take place, we are interested in the line 
density of zeros per unit length of this curve. We introduce the arclength s which measures the 
distance along C from the transition point zq. 
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FIG. 1: Small area A covering a length 5s of the curve of zeros C required to derive a relationship for the 
line density /i(s) of zeros along the curve. 

To obtain an expression for the line density of zeros /x(s) along the curve C, consider a short 
section of C of length Ss enclosed by an area A that has two sides parallel to C and the other two 
sides perpendicular to C — see figure 1. Integrating the density p(z) over this area we have 

I ^zp{z) = ii{s)bs . (11) 
JA 

Meanwhile, we have from (9) and an application of the divergence theorem that 

\^^zp{z) = i- ^dzV • [V^(z)] = ^[V¥'2(^) - V<^i(z)] ■ n (12) 

in which the functions (^1(2;) and (^2(2) relate to the limiting values of (/?(z) as the point z on the 
curve C is approached from either side, and the vector n is the unit vector normal to C at that 
point. 

Recall that (^(2;) is the real part of the complex free energy h{z) and therefore, away from any 
zeros of the partition function, satisfies the Cauchy-Riemann equations 

dx^{z) = dyipiyz) and dyip{z) = -d^'^^z) (13) 

in which il}{z) is the imaginary part of h{z) and z = x + iy. Then we have that 

V(^(z) • h = Vip{z) ■ i (14) 

where t is the unit vector tangent to the curve C at the point z. Thus we recognise the scalar 
product in (12) as the directional derivative of the imaginary part of the free energy h{z) along 
the curve C. Putting this together with (11) we find that 

l^{s) = :i^[Mz)-Mz)] (15) 
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in which iIj{z) = lm.[h{z)] and the subscripts indicate opposite sides of a phase boundary. 

Let us now assume that in the thermodynamic hmit, there is a phase transition at the point zq. 
Then, either side of the transition the free energy can be written as 

fl,2{z) = h{zQ) + bi,2Z + Ci,25^ + . . . (16) 

where z = z — zq and fi{z) is vahd for Rez < and f2{z) valid for Rez > 0. Note that for the free 
energy to be real along the real z axis, the coefficients b and c must also be real. Prom the criterion 
that the real part of h{z) is continuous across a phase boundary, we find that the boundary lies 
along the curve 

C2 - Ci 

where x and y are the real and imaginary parts of z respectively. We consider the conditions under 
which a discontinuous (first order) or continuous (second or higher order) appear. 

First-order transition 

For the case where bi ^ b2, and the free energy has a discontinuity in its first derivative, the curve 
of zeros is a hyperbola that passes smoothly through the transition point zq. Hence the tangent to 
the curve of zeros is parallel to the imaginary axis at zq and using the rule (15) we find that 

62 - 61 



Ai(0) = ^ ^[(62 - bi)y + 2(c2 - ci)xy] 



x=y=o 27r 



(18) 



Hence we see that the density of zeros at the transition point zq is nonzero at a first-order phase 
transition (i.e. one at which the first derivative of the free energy is discontinuous). 

Second-order transition 

If bi = 62 but C2 7^ ci we have that the curve C obeys the equation y = ±x. Since the zeros of the 
partition function come in complex-conjugate pairs, we find that the zeros approach the point zq 
along straight lines that meet at a right-angle. If we consider the line x = y = s/\/2, we find from 
(15) that 

This reveals that at a second-order phase transition, the density of partition function zeros decreases 
linearly to zero at the phase transition point. 

Higher-order transition 

More generally, one can consider a difference in the free energies either side of the transition point 
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to have the leading behaviour /2(n) — fi{u) ~ n". Then the condition that Re[f2{u) — fi{u)] = 
suggests that the curve of zeros approaches the real axis at an angle Unless a = 1, the curve 
does not pass smoothly through the real axis. In any case, the imaginary part of the free energy 
difference grows as giving rise to a density of zeros that behaves as for small arclength s. 
This means that for the density of zeros to be finite at the transition point we must have a > 1. 

In this section we have summarised what we refer to as the Lee- Yang theory of phase transitions. 
This describes how partition function zeros are related to phase transitions: the accumulation of 
zeros at a point along the physical (real, positive) axis of the control parameter gives the critical 
value and the density of zeros near to the accumulation point determines the order of the phase 
transition. 

We derived a rule (10) for locating phase boundaries and a further rule (15) for finding the 
density of zeros along such boundaries. At a first-order transition there is a nonzero density of 
zeros at the transition point whereas the density decays as a power-law to zero at the transition 
point when the associated phase transition is continuous. 

Although we have discussed the theory of partition function zeros with reference to the specific 
system described by the partition function (1) we should note that the ideas hold much more 
generally. Firstly, one is not restricted to the canonical ensemble: indeed in the original exposition 
of the theory [2], Lee and Yang worked in the grand canonical ensemble in which the quantity 
generalised to the complex plane was not a function of the temperature but rather the chemical 
potential. Of course, in the equilibrium theory, these intensive field-like quantities play similar 
mathematical roles and so there is no reason why the Lee- Yang theory shouldn't apply to all 
of them. For historical reasons, zeros in the complex fugacity (or chemical potential) plane are 
often referred to as Lee- Yang zeros [2, 3] and those in the complex temperature plane Fisher zeros 
[4]. Also it appears that one can just as easily generalise physical 'field-like' variables (such as 
temperature) or 'fugacity-like' variables (such as the quantity z considered above) to the complex 
plane without altering the properties of the partition-function zero densities at first-order and 
continuous phase transitions described above. 

Despite the apparently wide generality of the Lee- Yang theory of equilibrium phase transitions, 
proving its validity in the general case is a difficult task (although we note recent work in this 
direction [10]). Therefore, most rigorous results, such as those discussed in [8], tend to rely on 
specific properties of a particular partition function. Perhaps the most spectacular of these is that 
obtained by Lee and Yang in their original work. Specifically, they found that the zeros of the 
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partition function for an Ising ferromagnet in an external field h all lie on the circle | exp(/i)| = 1. 
The significance of this result, which does not depend on the number of spatial dimensions, lattice 
structure and details of the spin-spin interactions, is that if there is a phase transition induced by 
varying the magnetic field h, it can only occur at the point h = 0, and then only if the partition 
function zeros accumulate there in the thermodynamic limit. 

III. NONEQUILIBRIUM STEADY STATES 

A nonequilibrium steady state differs from its equilibrium counterpart in that it may admit 
a flow of, say, energy or mass. More generally, these states have a circulation of probability 
in the space of microscopic configurations. In the past few years, it has become customary to 
model nonequilibrium systems as stochastic processes, i.e. simple models defined by local dynamical 
rules. Extensive study of these has revealed a range of phenomena including nonequilibrium phase 
transitions [11-17]. We present below the key elements of this approach to nonequilibrium physics 
with a view to understanding nonequilibrium phase transitions in the framework of the Lee- Yang 
theory described above. In particular, we will need to propose a quantity to use in place of the 
equilibrium partition function (1). 

Let us begin by discussing how one might realise the dynamics of the equilibrium spin system 
of the previous section, for example in a computer simulation. The aim is to generate a sequence 
of spin configurations such that the frequency with which a particular configuration C is gener- 
ated is proportional to the Boltzmann weight f{C) = exp[— /3£'(C)] where E{C) is the energy of 
configuration C. Usually this is achieved by choosing the next configuration C' from the current 
configuration C with a probability proportional to the transition rate W(C C) that satisfies the 
detailed balance condition 

f{c)wic ^C) = fic')wic' ^c) yc,c' ■ (20) 

Apart from the fact that this choice of transition rates guarantees convergence to the desired 
equilibrium distribution of microstates [18], it also has the feature that the mean current of any 
quantity in the steady state is zero (as it should be at equilibrium). 

Since we are interested in nonequilibrium systems that support steady-state currents we must 
work with transition rates that do not satisfy the condition (20) or equivalent criteria (given in 
e.g. [14, 17, 19]). Clearly one has a lot of freedom in the choice of transition rates, and so in 
practice one often devises dynamics that seem physically reasonable given the phenomena one is 
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trying to describe. Later, in section IV A we will give concrete examples in the form of driven 
diffusive systems, in which particle moves are biased in a particular direction to model the effect 
of an external field. 

In order to find the set of steady-state weights associated with a particular set of prescribed 
transition rates, wc impose the condition that the total flow of probability into a configuration C 
is balanced by the corresponding outflow. That is, we must have 

[f{C')W{C' ^ C) - f{C)W{C ^ C')] = (21) 

C'^C 

for every configuration C. In general there might be more than one solution to this set of equations 
for the rates /(C), each corresponding to a steady state that is reached with some probability that 
depends on the starting configuration. We shall assume for simplicity that the steady state is 
unique and therefore reached with certainty from every initial state. A sufficient condition for a 
unique steady state is that it is possible to reach each configuration from every other via a sequence 
of microscopic transitions [17]. 

Once the steady state weights f{C) are known from (21) one can obtain the corresponding 
probability distribution of microscopic configurations P{C) through a normalisation 

Z = Yf{C) (22) 
c 

such that P{C) = f{C)/Z. Then, one can compute averages of physical quantities and look for 
nonanalyticities to locate nonequilibrium phase transitions as one varies the transition rates. At 
equilibrium we saw that the zeros of the function (1) that normalises the Boltzmann weights 
encodes the phase behaviour of the system. Thus we might hope that more generally, the zeros of 
the normalisation (22) will provide information about nonequilibrium phase behaviour. 

In section IV we review recent work that suggests that the steady-state phase behaviour of 
certain nonequilibrium driven diffusive and reaction-diffusion systems is correctly described by the 
zeros of the normalisation Z. Although we are unaware of any rigorous argument for this to be 
the case, some suggestive evidence is provided by an observation made in [17, 20] which we now 
outline. 

First note that the equation (21) is linear in the weights f{C). This implies that these weights, 
and hence the normalisation, can always be written as sums of products of the transition rates 
W{C C). Now, by numbering the microscopic configurations, one can construct the transition 
rate matrix W whose off-diagonal elements Wnm are equal to the transition rates from configuration 
m to n and the diagonal elements W„n are negative and express the total rate of departure from 
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configuration n. Using elementary results from matrix theory (see [17, 20] for the details), one 
can show that an expression for the normalisation Z that is polynomial in the transition rates 
W(C ^ C) is 

Z=ll{-\i) (23) 

in which Aj are the eigenvalues of the matrix W. 

With each of the eigenvalues Aj is associated an eigenvector of W describing a 'mode' of the 
stochastic process that decays exponentially with a timescale = l/|Ai| [17]. We have assumed 
that the process described by W has only one steady state, and so only one of the eigenvalues is 
equal to zero (since clearly the relaxation time of a steady state is infinite). Equation (23) states 
that the normalisation Z can be written as a product of the remaining eigenvalues. Now, at a phase 
transition we expect diverging timescales: physically one encounters metastable (long-lived) states 
near first-order transition points and long correlation lengths and times at continuous transitions. 
The presence of long timescales implies small eigenvalues Aj of W and, from (23), that the 
normalisation Z approaches zero. Hence it appears that it is appropriate to consider the zeros of 
Z as given by (22) to locate nonequilibrium phase transitions. 

IV. APPLICATION OF LEE- YANG THEORY TO NONEQUILIBRIUM PHASE 

TRANSITIONS 

We shall now review recent progress in applying the Lee- Yang theory to nonequilibrium phase 
transitions. We consider first in section IV A driven diffusive systems, where most work has so 
far been focussed [21-23]. An appealing feature, discussed in more detail below, is that some one- 
dimensional cases have been solved exactly, the normalisation (23) calculated, and nonequilibrium 
phase transitions analysed. (The existence of one-dimensional phase transitions contrasts with the 
case of one-dimensional equilibrium models which do not admit phase transitions if the interactions 
are short-ranged.) Then in Section IV B we move onto reaction-diffusion systems and directed 
percolation. 

A. Driven diffusive systems 



In their original papers on partition function zeros, Lee and Yang [2, 3] made use of the mapping 
between the Ising model of a magnet and the lattice gas. Essentially one associates up-spins with 
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particles and down-spins with vacancies so that a positive interaction strength J in the Ising model 
results in a particle-particle attraction in the lattice gas, and a negative interaction strength gives 
rise to repulsive gas particles. Note that there is implicitly at most one particle per lattice site, so 
there is a hard-core exclusion in the lattice gas. 

As discussed in section III, one can realise the dynamics of the lattice gas through a set of 
transition rates that satisfy the detailed balance condition (20) with respect to the Boltzmann 
distribution. Thirty years after Lee and Yang's work, Katz, Lebowitz and Spohn (KLS) [24, 25] 
introduced a driven lattice gas model in which the rate at which particles hop in the direction 
of an external field is enhanced and the hop-rate against the field is suppressed. This model is 
well-studied and many results are discussed in [26]. The principal effect of the driving field is to 
introduce anisotropy in the phase-separation that occurs below the critical temperature associated 
with the spontaneous magnetisation of the Ising ferromagnet (J > 0) in two or more spatial 
dimensions and in the critical exponents that characterise this transition. 

As yet, the KLS model remains unsolved for general interaction strength, although in one 
dimension the steady state is known for some parameters [27] . The particular case of one dimension 
and zero interaction strength J = 0, known as the asymmetric simple exclusion process (ASEP), 
had already been studied — at least at a mean-field level — by biophysicists interested in the kinetics 
of biopolymerisation [28] . The mean-field approach predicts phase transitions in the steady state 
as parameters controlling the rate of insertion and extraction of particles at the boundaries are 
varied [29]. The existence of these phase transitions is confirmed through an exact solution of the 
ASEP [29-31], achieved using a powerful matrix product approach [30, 32] which has subsequently 
been used to solve many generalisations of the ASEP. The details of the matrix product method 
are not necessary for the following — suffice to say that one ends up calculating a normalisation 
proportional to (23) through a product of matrices, often of infinite dimension. In this way one 
obtains some explicit formulas for the normalisation (22) which we shall use below [50] . 

The asymmetric exclusion process with open boundaries is perhaps the simplest exactly solved 
nonequilibrium model that exhibits both a first-order and continuous phase transition in its steady 
state. Therefore it is an ideal candidate for testing the hypothesis outlined in section III that 
zeros of the normalisation should accumulate towards the positive real axis in the complex plane 
of transition rates. Before outlining the results of this analysis (the details of which are presented 
in [22]) we recall the definition of the ASEP with open boundaries. 

In this system, a particle on an A?^-site lattice can hop one site to the right at unit rate, as long 
as the receiving site is empty. Meanwhile, particles are inserted onto the leftmost lattice site (if 
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FIG. 2: The dynamics defining the ASEP. The lattice is open and particles may be inserted or extracted at 
the boundaries as shown. 
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FIG. 3: The phase diagram of the ASEP in the space spanned by the boundary rates a and /3. The thick 
line represents a first-order phase boundary, the thin line a continuous phase boundary. 

empty) at a rate a and removed from the rightmost site (if occupied) at a rate /? — see figure 2. 
Along the line a = (3 < ^ there is a coexistence between a high- and a low-density phase, with 
a shock separating the two. This is indicative of a first-order transtion. Meanwhile, along the 
lines a > = ^ and (3 > = ^ there is a continuous transition (i.e. one accompanied by 
diverging lengthscales) to a phase in which the particle current is independent of a and (3. The 
phase diagram for the model is shown in figure 3. 

The steady-state normalisation (22) for the ASEP with N sites has been calculated [30] as 

^"(".ffl^g jvKiv-rt! - (1/.) ■ 

It is a simple matter to use a computer algebra package to solve this equation for its zeros in 
the complex-a plane at fixed N and /3. (Equivalently, one could look at the complex-/? zeros at 
fixed a since Z]sf{a,l3) = Zi^{f3,a).) In figure 4 plots of the zeros are shown for /? = 1, where 
a continuous transition occurs at a = ^, and for /3 = |, where a first order transition occurs at 
a = \. We immediately notice that the curve of zeros seem to intersect the real positive a axis at 
the correct transition point. Furthermore, the density of zeros near the first-order transition point 
{a = (3 = \) seems to be uniform and nonzero, whereas the density of zeros near the second-order 
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FIG. 4: Zeros of the normalisation in the ASEP in the complex-a plane and /3 = ^,1 and lattice size 
iV = 300. 

transition point (a = ^, /? = 1) seems to decrease to zero. Both of these observations are in accord 
with the results known for equilibrium partition function zeros discussed in section II. 

As shown in [22] , the distribution of zeros in the thermodynamic limit can be calculated once 
one knows that, for large AT, the normalisation behaves as Zj^ ~ AJ-^N'i. In this expression, 
A, J and 7 depend on a and /3 and the quantity J is the current of particles across the lattice. 
Thus, the complex free energy (3) in the limit AT — oo is 

= lim ^^^^("'/^) = - In J(a,/3) . (25) 

AT— MX) A* 

Furthermore for a smaller than the transition value ac-, J = cai}- — ol) whereas ioY a > ac-, J = 
ac(l - Oc)- 

Although an electostatic analogy was used in [22] to find the zero distribution, the mathematical 
content is the same as that used to derive the two rules (10) and (15) in section II. Applying the first 
rule, which demands that the real part of the free energy be continuous across a phase boundary, 
we find that the zeros of Z^{a^ /?) should lie along the curve |a(l — a)| = aJ^V — ac). It is therefore 
convenient to change variable to ^ = a(l — a). Then, in the complex-^ plane, the zeros lie on the 
circle |^| = = ctc(l — etc)- The density of zeros on this circle can be found by setting ^ = re*^ 
and parametrising points on the circle as s = ^c^- Then, the second rule (15) gives the density of 
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zeros fx{s) on the circle as 

That is, in the complex-,^ plane the zeros should become evenly distributed on a circle in the 
thermodynamic limit. Transforming the zeros of (24) obtained at different system sizes to the 
complex-^ plane reveals this to be the case [22]. 

Finally, one can show that near the intersection point between the curve of zeros in the complex- 
a plane and the positive real a axis, the zeros of (24) sit on the curve x = ^ — {y'^ + \ — ^c)^^^ 
where x and y are the real and imaginary parts of a. For the case /? < 5, < j and the transition 
is first order. One finds that the curve of zeros is smooth at the transition point a = (3, and that 
the density of zeros is (1 ^ 2/3)/[27r/3(l — (3)] which is nonzero. These are precisely the properties 
of the equilibrium partition function zeros at a first-order transition point (see section II). 

At the continuous transition point (/? > ^ and £, = \, the zeros pass through the real a axis 
along the line x = ^ — \y\. Recall from section II that a transition of n}^ order has the curve of 
zeros meeting the positive real axis at an angle of ^ . Here the zeros clearly approach at an angle 
I suggesting a second-order transition. In fact, one finds the density of zeros is |s at a distance s 
along this line, confirming that the transition is second-order. 

In summary, we have found that the Lee- Yang theory of first-order and continuous phase tran- 
sitions applies to the normalisation of the nonequilibrium asymmetric exclusion process just as 
it does to the partition function of equilibrium systems. Of course, this does not prove that the 
theory is generally applicable, and so there is some value in investigating other nonequilibrium 
steady states that exhibit phase transitions. 

One such state is that initially studied by Arndt, Heinzel and Rittenberg [33, 34]. This model 
comprises M+ (M_) positively (negatively) 'charged' particles on a closed ring of N sites. When 
next to vacant sites, the positive particles hop to the right and negative particles to the left at unit 
rate, with hops in the opposite directions disallowed. Meanwhile, should two oppositely charged 
particles be next to one another, the following transitions can occur 

where the label above the arrow indicates the rate at which the transitions occur. These dynamics 
are illustrated in figure 5. 

Although a matrix product solution for this model is known to exist [34] , it is technically difficult 
to work out the solution in its full generality [35] . However in the case of either a single vacancy 




FIG. 5: The dynamics defining the AHR model: labels indicate the rates at which the hops denoted by 
arrows can occur. Note that the N-site lattice is periodic. 

[36] or a single negative particle, i.e. M_ = 1 [37] the steady-state normalisation is known explcitly 
and phase transitions have been identified. In the case M_ = 1 the zeros of the partition function 
in the complex-g plane have been studied [23] and, once again the continuous phase transition 
present in this model is correctly described by the Lee- Yang theory. 

A slightly different Lee- Yang analysis, that precedes the work described above, has been used 
to study the nonequilbrium phase transition that occurs at a finite density of vacancies [21]. The 
motivation behind this study arose from computer simulations of the model [33, 34] for the case 
where the numbers of positive and negative particles were equal. In these simulations, three phases 
were initially identified: a 'solid' phase for small q, a 'fluid' phase for large q and a mixed phase 
comprising the background 'fluid' and a large mobile droplet in an intermediate regime 1 < q < qc 
where qc is some density-dependent quantity. 

The steady-state normalisation Zjv,m for the case of = M_ = M positive and negative 
particles on the iV-site ring is not known exactly (at least in the canonical ensemble), therefore 
Arndt [21] chose to study the generating function 

ZN{z,q)=J2^^^^M<l) (27) 
M 

in the complex-z plane. Physically this approach is equivalent to placing the ring with nonequi- 
librium interactions in chemical equilibrium with a particle reservoir at fixed fugacity z. One can 
check that as one takes the thermodynamic limit, the relative size of the fluctuations in the density 
r = M/N vanish, and so working at fixed fugacity z is equivalent to working at fixed density r. 
(In fact, this is a standard 'trick' for dealing with closed systems in which the particle number is 
conserved, sec e.g. [38].) 

The key point here, however, is that z is an equilibrium fugacity, and not a microscopic transition 
rate, and so the Lee- Yang theory of phase transitions described in section II ought to apply directly 
here, without reference to the discussion of section III. It was found [21] that when q < qc the 
zeros of (27) in the complex-z plane appear to lie along ellipses, intersecting the positive real axis. 
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Conversely, for larger q the zeros appear to describe hyperbolae that avoid the positive real axis. 
After investigating numerically the density of zeros in the complex-z plane, Arndt concluded that 
there is a fluid-fluid phase transition at some density-dependent point > 1- 

Unfortunately, an exact asymptotic (i.e. large N) analysis of the generating function (27) shows 
that the only nonanalyticity occurs at the solid-fluid transition point q = 1 [35]. However it 
was noted that physical quantities vary rapidly, but continuously, at some a point g > 1. This 
phenomenom has been explained as an abrupt increase in a corrleation length to an anomalously 
large, but finite value [39]. We believe that resolution with the study of Arndt lies in the fact 
that the system sizes considered in [21] were quite small {N < 100), whereas it was suggested [35] 
that the distinction between this crossover behaviour and a genuine phase transition might become 
apparant only once N is increased above about 10^° (an unfeasibly huge number!). 

It would be interesting, therefore, to extend the numerical computation of zeros performed in 
[21] to much larger systems, to see how the ellipses noted above develop. For example it might 
well happen that instead of approaching the positive-real fugacity axis, the zeros of (27) would 
terminate a short distance away from it. 

B. Reaction-diffusion systems and directed percolation 

A large and important class of models with stochastic dynamics is provided by reaction-diffusion 
systems [51]. In contrast to driven diffusive systems, where the particle-particle interactions imply 
conservation of particles, reaction-diffusion systems are characterised by dynamics that result in 
a change in particle number. Moreover there are a number of such systems that have absorbing 
states i.e. special configurations generally devoid of particles that once entered cannot be left. Phase 
transitions associated with whether the system has a finite probability of not being absorbed into 
such a state fall within the directed percolation and related universality classes [13]. We shall shortly 
discuss the second order phase transition associated with the directed percolation university class 
in a little detail. Meanwhile as a simple example of a reaction-diffusion system, we review a model 
for which the steady state that can be solved using the matrix product approach. The approach 
again provides us with an explicit expression for the normalisation (22) by virtue of which we can 
analyse its zeros in the plane of complex reaction rates. 

The system in question [40] has for the dynamics at neighbouring bulk sites on a one-dimensional 
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lattice the processes 
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occurring at the rates indicated, and with • representing a particle and o an empty lattice site. It 
was demonstrated [40] that the matrix product scheme used for the ASEP could be generalised 
to cater for the steady state of the present reaction-diffusion system on a lattice of N sites with 
reflecting boundary conditions (i.e. which is neither periodic nor has particle input or removal at 
the boundaries). 

In the matrix product approach, the normalisation is given by a scalar derived from 
where C is a square matrix. If C is a finite dimensional m x m matrix and can be diagonalised, 
the resulting expression for has the form 

m 

ZN{q,n)=Y,anXn (28) 
n=l 

in which A„ is the n*^ eigenvalue of C and the coefficients a„ arise from the details of the way in 
which one obtains a scalar from the matrix . A normalisation of this form leads to a complex 
free energy 

h{q, k) = In max{An} (29) 

In J 

in which the maximum means the eigenvalue with the largest absolute value. Then, there is the 
possibility of a phase boundary when the magnitude of the two largest eigenvalues of C are equal. 
It turns out that for the closed reaction-diffusion system introduced above, C is a 4 x 4 matrix that 
can be diagonalised when ^ 1 + k [40]. When q'^ = 1 + k, C cannot be diagonalised on account 
of the largest eigenvalue being degenerate, and a phase transition is found to occur at this point. 
Note that this scenario contrasts with the transfer-matrix approach to one-dimensional equilibrium 
systems where the partition function is also written as a product of matrices. Since all elements 
of the transfer matrix arc positive the largest eigenvalue cannot become degenerate therefore there 
can be no phase transition [12]. However there is no such restriction on the elements of C, and so 
eigenvalue crossing is permitted and nonequilibrium one-dimensional phase transitions can occur. 

The structure of the normalisation Zjv(q, k) is even simpler in the case where particles are 
inserted onto and removed from the left boundary at rates a and (3 respectively (with the right 
boundary remaining reflecting), and these rates satisfy the relation a = K{q~^ — q + P) [41]. In this 
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FIG. 6: Directed percolation. Bonds are open (solid lines) with probability p and closed (broken lines) with 
probability 1 — p. A fluid starts at layer n = and may flow only downwards through open bonds. Bonds 
carrying the fluid are shown with thick lines, so that the point A on layer n = 6 is connected to the origin 
whereas point B is not. 

case C is a 2 X 2 matrix and, for 7^ 1 + k, is diagonalisable with eigenvalues 



We see then immediately from (29) and an application of the rule (10) that there is a phase 
boundary in the complex q plane along the circle l^l = \/l + 1^- Also, since one of the eigenvalues 
does not depend on the free energy /i(g, k) is a constant in one of the phases which, as the 
analysis of the normalisation zeros for the ASEP demonstrated, implies that the density of zeros 
on this circle is constant in the thermodynamic limit. In turn, this implies that the phase transition 
is first order, as confirmed by explicitly calculating the density profile in the two phases [41]. 

We finally turn our attention to a process that comprises symmetric decoagulation (that is, 
•o — > and o» •• taking place with equal probability in each direction) and spontaneous decay 
(• — o) occurring at independent rates. Introduced as a crude model of an epidemic [42], this 
contact process is known to exhibit a transition from a phase in which the absorbing state (empty 
lattice) is reached with certainty to a phase in which there is some probability that the epidemic 
remains active forever in the thermodynamic (infinite system size) limit [43] . This transition occurs 
as the ratio between the decoagulation and decay rates is increased beyond a critical value. 

Although not proven, it is widely accepted on the basis of simulation and approximate methods 
(such as series expansions) that the transition just described is continuous and characterised by 
directed percolation (DP) exponents (see [13] for an in-depth discussion of these issues). Directed 
percolation itself [44] is a geometric construction designed to model fluid flow through a random 
porous medium under the influence of gravity. Consider the rhombic lattice shown in figure 6 



Ai = 1 + At and X2 = ■ 



(30) 
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and imagine 'pouring' fluid in from the uppermost site. The fluid is allowed to flow along a bond 
pointing diagonally downwards and then only if it is open. Each bond has a probability p of being 
open (and hence a probability 1 — p of being blocked), and the state of each bond is independent 
of any other and may not change with time. 

The main quantity of interest in this system is the percolation probability Pn{p) that the 
fluid can penetrate to a depth n in the medium (see figure 6), averaged over all possible bond 
configurations. The order parameter for the model is the probability Pooip) that the fiuid can 
penetrate infinitely far given a bond probability p. If p is below some critical percolation threshold 
Pc, the fluid only ever penetrates a finite distance — i.e. a layer becomes dry with certainty. However, 
for larger p, the order parameter becomes nonzero, growing as Poo{p) ~ \p — Pc\^ with P 0.276 
near the critical point. The importance of the DP transition is that a wide range of models that 
have a transition into an absorbing state are expected to have that transition characterised by the 
DP exponents, one of which is /3 [13]. Despite a huge amount of interest in directed percolation, 
none of the models expected to belong to its universality class has been solved exactly. 

Recently an attempt has been made to shed further light on the DP transition by studying 
the zeros of the percolation probability Pn{p) [45, 46]. At first glance, a connection between this 
probability and a partition function is not obvious. However, associating with each site a 'spin' 
state CTj, that is up if site i is connected to a point on the n*^ layer and down otherwise, yields a 
form for Pn{p) that has the structure of a transfer-matrix representation of a partition function for 
an equilibrium system with three-spin interactions [47]. 

It is not clear that being able to give Pn{p) the appearance of a partition function necessarily 
implies that the Lee- Yang theory should hold. In particular, Pn{p) has some features that set it 
apart from 'standard' equilibrium partition functions, in that it does not grow exponentially in 
the number of bonds at least in the range < p < 1. It is also uncertain whether a free 
energy can meaningfully be defined for this system, especially in the region < p < Pc in which 
limn^oo Pn{p) = 0. 

Given these observations, it is perhaps not surprising that the zeros of Pnip) calculated nu- 
merically for n < 15 (leading to polynomials of degree up to = 240) have a more complicated 
distribution than for the (exactly solved) models considered so far. As is evident from figure 7, 
which was presented in [46] , the zeros do not lie along a single curve but along a sequence of curves 
meeting at the critical point (as well as inside some region that encloses part of the real axis for 
P>'2). 

By considering the sequences of zeros approaching the positive real axis, one can estimate the 
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FIG. 7: Zeros of the percolation probability on the directed percolation lattice up to depths n = 15. Taken 
from [46] in which the symbol t is used where we use n. 

critical point Pc and the density of zeros as it is reached [46] . The latter yields a prediction for one 
of the DP critical exponents and one finds good agreement with the most precise estimates of both 
the transition point and the associated exponent obtained through other means. Hence we have 
further evidence for the applicability of Lee and Yang's ideas concerning partition function zeros 
to a much wider range of statistical distributions than equilibrium steady states. 



V. SUMMARY AND OUTLOOK 



In this work we have revisited the Lee- Yang description of equilbrium phase transitions with a 
view to seeing whether the ideas apply to more general nonequilibrium transitions. Recently there 
have been a number of studies of zeros of partition-function-like quantities that arise in systems 
with nonequilibrium dynamics, and we have seen in our review of these works that the Lee- Yang 
theory, as described in section II, seems to hold quite generally. 

We have argued that for dynamic models with a unique steady state, the normalisation defined 
as a sum over the steady-state configurational weights (22) serves as a suitable 'partition function' 
in the sense that its zeros, in the complex plane of any model parameter, should accumulate towards 
physical transition points in the thermodynamic limit. Furthermore, the density of zeros and angle 
of approach to the real axis indicate whether the transition is first-order (manifested physically 
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through phase-coexistence) or continuous (i.e. characterised by divergent correlation lengths and 
times). Thus studying the zeros of the normahsation (22) provides an unambiguous classification 
of nonequilibrium phase transitions as do the Lee- Yang zeros in the equilibrium case. 

The observation that backs up this scenario is embodied by equation (23) which reveals that the 
reciprocal of the steady-state normalisation is equal to the product of the characteristic relaxation 
times in the dynamics. Since near a phase transition one expects timescales to diverge, one also 
expects the normalisation to approach zero. However, a rigorous argument for this to be the case 
is still lacking. Moreover (23) implies a possible link between systems for which the steady state 
normalisation can be calculated and those for which eigenvalues of the transition matrix can in 
principle be calculated. 

A different class of systems encompasses those whose steady state is not unique. The contact 
process is, in fact, an example of such a model, in which the absorbing state is reached with certainty 
below the critical dccoagulation rate, whereas above it, and on an infinite system, a second steady 
state can also be reached with some nonzero probability. Since this additional steady state exists 
only when the lattice size becomes infinite, one must take that limit first, before taking time to 
infinity. Otherwise, on a finite system the steady state is simply the absorbing state and the 
steady-state normalisation is trivially equal to a constant. Nevertheless the work of [45, 46], which 
we reviewed in section IV B, indicates that the Lee- Yang theory can be relevant when one considers 
other properties of the nonequilibrium system such as the percolation probability. 

Although we have not discussed this in great detail here, it should be noted that the Lee- Yang 
approach gives a method for extrapolating to the thermodynamic limit from solutions for small 
system sizes. In the work of [46], mimerical solutions for small systems were used successfully to 
estimate the transition point and density of zeros as it is approached. Prom this information one 
learns about the nature of the phase transition and, for example, can estimate the values of critical 
exponents. This is a common technique in equilibrium statistical physics (see e.g. [48]) and it may 
be the case that stochastic processes unyielding to analytical treatment could be understood this 
way. 



22 



Acknowledgments 

We thank Bernard Derrida for helpful discussions. R.A.B. acknowledges financial support under 
EPSRC grant GR/R53197. 



[1] C. Domb (1996) The critical point: a historical introduction to the modern theory of critical phenomena 

(Taylor & Francis, London). 
[2] C. N. Yang and T. D. Lcc (1952) Statistical theory of equations of state and phase transitions. I. Theory 

of condensation. Phys. Rev. 87 404. 
[3] T. D. Lee and C. N. Yang (1952) Statistical theory of equations of state and phase transitions. IL 

Lattice gas and Ising model. Phys. Rev. 87 410. 
[4] M. E. Fisher (1965) The nature of critical points. In Lectures in theoretical physics, edited by W. E. 

Brittin, volume 7C, page 1 (University of Colorado Press, Boulder). 
[5] S. Grossmann and W. Roscnhauer (1969) Phase transitions and the distribution of temperature zeros 

of the partition function: I. General relations. Z. Physik 218 437. 
[6] S. Grossmann and V. Lchmann (1969) Phase transitions and the distribution of temperature zeros of 

the partition hmction: II. Applications and examples. Z. Physik 218 449. 
[7] L. E. Reichl (1980) A modern course in statistical physics (University of Texas Press, Austin). 
[8] D. Ruelle (1969) Statistical mechanics: Rigorous results (W A Benjamin, New York). 
[9] O. D. Kellogg (1953) Foundations of Potential Theory (Dover, New York). 
[10] M. Biskup, C. Borgs, J. T. Chaycs, L. J. Klcinwaks and R. Kotccky (2000) A general theory of Lee- Yang 

zeros in models with first-order phase transitions. Phys. Rev. Lett. 84 4794. 
[11] .1. Marro and R. Dickman (1999) Nonequilihrium, phase transitions in lattice models (CUP, Cambridge). 
[12] M. R. Evans (2000) Phase transitions in one-dimensional nonequilibrium systems. Braz. J. Phys. 30 

42. 

[13] H. Hinrichscn (2000) Nonequilibrium critical phenomena and phase transitions into absorbing states. 
Adv. Phys. 49 815. 

[14] D. Mukamel (2000) Phase transitions in nonequilibrium systems. In Soft and fragile matter, edited by 

M. E. Gates and M. R. Evans, volume 53 of Scottish Universities Summer School in Physics, page 237 

(Institute of Physics Publishing, Bristol). 
[15] G. M. Schiitz (2001) Exactly solvable models for many-body systems far from equilibrium. In Phase 

transitions and critical phenomena, edited by C. Domb and J. L. Lebowitz, volume 19, page 1 (Academic 

Press, London). 

[16] R. Stinchcombe (2001) Stochastic non-equilibrium systems. Adv. Phys. 50 431. 

[17] M. R. Evans and R. A. Blythe (2002) Nonequilibrium dynamics in low-dimensional systems. Physica 



23 



A 313 110. 

O. Narayan and A. P. Young (2001) Convergence of Monte-Carlo simulations to equilibrium. Phys. 
Rev. E 64 021104. 

N. G. van Kampen (1992) Stochastic processes in physics and chemistry (North-Holland, Amsterdam), 
second edition. 

R. A. Blythe (2001) Nonequilibrium phase transitions and dynamical scaling regimes. Ph.D. thesis, 
University of Edinburgh, http://www.ph.ed.ac.uk/cmatter/links/rab-thesis. 

P. F. Arndt (2000) Yang-Lee theory for a nonequilibrium phase transition. Phys. Rev. Lett. 84 814. 
R. A. Blythe and M. R. Evans (2002) Lee- Yang zeros and phase transitions in nonequilibrium steady 
states. Phys. Rev. Lett. 89 080601. 

F. H. Jafarpour. The application of the Yang-Lee theory to study a phase transition in a non-equilibrium 
system. Preprint, cond-mat/0209558. 

S. Katz, J. L. Lebowitz and H. Spohn (1983) Phase transitions in stationary nonequilibrium states of 
model lattice systems. Phys. Rev. B 28 1655. 

S. Katz, J. L. Lebowitz and H. Spohn (1984) Nonequilibrium steady-states of stochastic lattice gas 
models of fast ionic conductors. J. Stat. Phys 34 497. 

B. Schmittmann and R. K. P. Zia (1995) Statistical mechanics of driven diffusive systems, volume 17 
of Phase transitions and critical phenomena (Academic Press, London). 

J. Krug and H. Spohn (1992) Kinetic roughening of growing surfaces. In solids far from equilibrium, 
edited by C. Godreche (Cambridge University Press, Cambridge). 

C. T. MacDonald, J. H. Gibbs and A. C. Pipkin (1968) Kinetics of biopolymerization on nucleic acid 
templates. Biopolymers 6 1. 

B. Derrida, E. Domany and D. Mukamel (1992) An exact solution of a one-dimensional asymmetric 
exclusion model with open boundaries. J. Stat. Phys. 69 667. 

B. Derrida, M. R. Evans, V. Hakim and V. Pasquier (1993) Exact solution of a ID asymmetric exclusion 
model using a matrix formulation. J. Phys. A: Math. Gen. 26 1493. 

G. Schiitz and E. Domany (1993) Phase transitions in an exactly soluble one-dimensional exclusion 
process. J. Stat. Phys. 72 277. 

B. Derrida (1998) An exactly soluble non-equilibrium system: The asymmetric simple exclusion process. 
Phys. Rep. 301 65. 

P. F. Arndt, T. Heinzel and V. Rittenberg (1998) Spontaneous breaking of translational invariance in 
one-dimensional stationary states on a ring. J. Phys. A: Math. Gen. 31 L45. 

P. F. Arndt, T. Heinzel and V. Rittenberg (1999) Spontaneous breaking of translational invariance and 
spatial condensation in stationary states on a ring: I. The neutral system. J. Stat. Phys. 97 1. 
N. Rajewsky, T. Sasamoto and E. R. Speer (2000) Spatial Particle Condensation for an Exclusion 
Process on a Ring. Physica A 279 123. 

T. Sasamoto (2000) One-dimensional partially asymmetric simple exclusion process on a ring with a 



24 



defect particle. Phys. Rev. E 61 4980. 
[37] F. H. Jafarpour (2000) Exact solution of an exclusion model in the presence of a moving impurity on 

a ring. J. Phys. A: Math. Gen. 33 8673. 
[38] B. Derrida, S. A. Janowsky, J. L. Lebowitz and E. R. Speer (1993) Exact solution of the totally 

asymmetric simple exclusion process: shock profiles. J. Stat. Phys. 73 813. 
[39] Y. Kafri, E. Levine, D. Mukamel and J. Torok (2002) Sharp crossover and anomalously large correlation 

length in driven systems. J. Phys. A: Math. Gen. 35 L459. 
[40] H. Hinrichsen, S. Sandow and I. Peschel (1996) On matrix product ground states for reaction - diffusion 

models. J. Phys. A: Math. Gen. 29 2643. 
[41] F. H. Jafarpour. First order phase transition in a reaction-diffusion model with open boundary: the 

Yang-Lee theory approach. Preprint, cond-mat/0301407. 
[42] T. E. Harris (1974) Contact interactions on the lattice. Ann. Prob. 2 969. 

[43] D. Griffeath (1979) Additive and cancellative interacting particle systems, volume 724 of Lecture notes 

in mathematics (Springer- Verlag, Berlin). 
[44] S. R. Broadbent and J. M. Hammersley (1957) Percolation processes I: Crystals and mazes. Proc. 

Camb. Phil. Soc. 53 629. 

[45] P. F. Arndt, S. R. Dahmen and H. Hinrichsen (2001) Directed percolation, fractal roots and the Lee- 
Yang theorem. Physica A 295 128. 

[46] S. M. Dammer, S. R. Dahmen and H. Hinrichsen (2002) Yang-Lee zeros for a nonequilibrium phase 
transition. J. Phys. A: Math. Gen. 35 4527. 

[47] R. J. Baxter and A. J. Guttmann (1988) Series expansion of the percolation probability for the directed 
square lattice. J. Phys. A: Math. Gen. 21 3193. 

[48] W. Janke and R. Kenna (2002) Phase transition strengths from the density of partition function zeroes. 
Nucl. Phys. B Proc. Supp. 106 905. See also hep-lat/01 12032. 

[49] S. Redner (1997) Scaling theories of diffusion-controlled and ballistically-controUed bimolecular reac- 
tions. In Nonequilibrium statistical mechanics in one dimension, edited by V. Privman, chapter 1 
(Cambridge University Press, Cambridge). 

[50] As yet, there is no comprehensive review of applications of the matrix-product method, although one 
should consult [13, 15, 16, 32] and references therein 

[51] As an area of ongoing research, new results are emerging continuously and there is no up-to-date review 
that contains them all. Nevertheless, one can consult [13, 49] for an introduction. 



